{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The Lyapunov spectrum of the Lorenz system"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we present a calculation of the Lyapunov spectrum of the Lorenz system, using `TaylorIntegration.jl`. Our calculation involves evaluating the 1st order variational equations for this system; the variational equations are evaluated by `TaylorIntegration.jl` via automatic differentiation techniques implemented in `TaylorSeries.jl`; this means that we don't have to calculate manually the variational equations, `TaylorIntegration.jl` does that for us! This saves time, as well as typos, which sometimes are hard to trace, especially for systems with a large number of dependent variables.\n",
    "\n",
    "The Lorenz system is the ODE defined as:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "    \\begin{split}\n",
    "    \\dot x_1 &= \\sigma(x_2-x_1) \\\\\n",
    "    \\dot x_2 &= x_1(\\rho-x_3)-x_2 \\\\\n",
    "    \\dot x_3 &= x_1x_2-\\beta x_3 \n",
    "    \\end{split}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "where $\\sigma$, $\\rho$ and $\\beta$ are constant parameters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The setup"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we write a Julia function which evaluates (in-place) the Lorenz system equations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "lorenz! (generic function with 1 method)"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Lorenz system ODE:\n",
    "function lorenz!(t, x, dx)\n",
    "    dx[1] = σ*(x[2]-x[1])\n",
    "    dx[2] = x[1]*(ρ-x[3])-x[2]\n",
    "    dx[3] = x[1]*x[2]-β*x[3]\n",
    "    nothing\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we define the parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "45.92"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Lorenz system parameters\n",
    "#we use the `const` prefix in order to help the compiler speed things up\n",
    "const σ = 16.0\n",
    "const β = 4.0\n",
    "const ρ = 45.92"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The initial conditions, the initial and final time are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "100.0"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "const x0 = [19.0, 20.0, 50.0] #the initial condition\n",
    "const t0 = 0.0 #the initial time\n",
    "const tmax = 100.0 #final time of integration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We know that the sum of the Lyapunov spectrum has to be equal to the trace of the Jacobian of the equations of motion. We will calculate this trace using `TaylorSeries.jl`, and after the numerical integration, we will come back to this value to check that this is indeed the case:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "using TaylorSeries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-21.0"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Calculate trace of Lorenz system Jacobian via TaylorSeries.jacobian:\n",
    "xi = set_variables(\"δ\", order=1, numvars=length(x0))\n",
    "x0TN = [ x0[1]+xi[1], x0[2]+xi[2], x0[3]+xi[3] ]\n",
    "dx0TN = similar(x0TN)\n",
    "lorenz!(t0, x0TN, dx0TN)\n",
    "lorenztr = trace(jacobian(dx0TN)) #trace of Lorenz system Jacobian matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The integration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we are ready to perform the integration, using `TaylorIntegration.jl`. Normally, we would use the `taylorinteg` method in order to integrate the equations of motion; but, since we are interested in evaluating the Lyapunov spectrum, we will use the `liap_taylorinteg` method, which calculates the Lyapunov spectrum via the variational equations and Oseledets' theorem. The expansion order will be $28$; the local absolute tolerance will be $10^{-20}$. `liap_taylorinteg` will return three arrays: one with the evaluation times, one with the value of the dependent variables at the time of evaluation, and another one with the vaues of the Lyapunov spectrum at each time of the evaluation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "using TaylorIntegration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " 11.295377 seconds (219.22 M allocations: 15.129 GiB, 15.68% gc time)\n"
     ]
    }
   ],
   "source": [
    "@time tv, xv, λv = liap_taylorinteg(lorenz!, x0, t0, tmax, 28, 1e-20; maxsteps=2000000);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The number of steps taken is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4218"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "length(tv)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Is the final time actually the requested value?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "true"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tv[end] == tmax"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What is the minimum, maximum and mean time-step?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.0018190222018432678, 0.04730949106561866, 0.023713540431586435)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "minimum(diff(tv)), maximum(diff(tv)), mean(diff(tv))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What is the standard deviation of the time-step size distribution?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.006662910504132029"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "std(diff(tv))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What is the final value for each of the exponents?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Float64,1}:\n",
       "   1.53183   \n",
       "  -0.00994121\n",
       " -22.5219    "
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "λv[end,:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What is the value of the sum of the value we obtained for the spectrum?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-20.99999999999998"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum(λv[end,:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, is the sum of the exponents exactly equal to the trace of the Jacobian?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "false"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum(λv[end,:]) == lorenztr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So the sum is not **exactly** equal, but is it approximately equal? What is the relative error in the computed vs the expected value?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-1.0150610510858574e-15"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rel_error = (sum(λv[end,:])-lorenztr)/lorenztr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How does this error compare to the machine epsilon?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-4.571428571428571"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rel_error/eps()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So the relative error in our computation for the sum of the Lyapunov spectrum is comparable to the machine epsilon."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Visualization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we plot our results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Plots.PyPlotBackend()"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using Plots\n",
    "pyplot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Lyapunov exponents vs time plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot(tv, λv[:,1], label=\"λ₁\")\n",
    "plot!(tv, λv[:,2], label=\"λ₂\")\n",
    "plot!(tv, λv[:,3], label=\"λ₃\")\n",
    "xlabel!(\"time\")\n",
    "ylabel!(\"λᵢ, i=1,2,3\")\n",
    "title!(\"Lyapunov exponents vs time\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Lyapunov exponents vs time plot (semi-log)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot(tv[2:end], λv[:,1][2:end], xscale=:log10, label=\"λ₁\")\n",
    "plot!(tv[2:end], λv[:,2][2:end], label=\"λ₂\")\n",
    "plot!(tv[2:end], λv[:,3][2:end], label=\"λ₃\")\n",
    "xlabel!(\"time\")\n",
    "ylabel!(\"λᵢ, i=1,2,3\")\n",
    "title!(\"Lyapunov exponents vs time (semi-log)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Convergence of Lyapunov exponents vs time plot (log-log)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot(tv[2:end], abs.(diff(λv[:,1])), yscale=:log10, xscale=:log10, label=\"δλ₁\")\n",
    "plot!(tv[2:end], abs.(diff(λv[:,2])), label=\"δλ₂\")\n",
    "plot!(tv[2:end], abs.(diff(λv[:,3])), label=\"δλ₃\")\n",
    "xlabel!(\"time\")\n",
    "ylabel!(\"δλᵢ, i=1,2,3\")\n",
    "title!(\"Convergence of Lyapunov exponents vs time (log-log)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Absolute difference wrt trace of Jacobian vs time plot (log-log)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot(tv[2:end], log10.(abs.(λv[2:end,1]+λv[2:end,2]+λv[2:end,3]-lorenztr)),label=\"log₁₀|Σᵢλᵢ-Tr(J)|\", xscale=:log10)\n",
    "ylims!(-15,-12)\n",
    "xlabel!(\"time\")\n",
    "ylabel!(\"log₁₀(absolute difference)\")\n",
    "title!(\"log₁₀|Σᵢλᵢ-Tr(J)| vs time (log-log)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Lorenz attractor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot(xv[:,1], xv[:,2], xv[:,3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.6.0-rc2",
   "language": "julia",
   "name": "julia-0.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
